Entanglement and density-functional theory: testing approximations on Hooke's atom. 
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We present two methods of calculating the spatial entanglement of an interacting electron system within 
the framework of density-functional theory. These methods are tested on the model system of Hooke's atom for 
which the spatial entanglement can be calculated exactly. We analyse how the strength of the confining potential 
affects the spatial entanglement and how accurately the methods that we introduced reproduce the exact trends. 
We also compare the results with the outcomes of standard first-order perturbation methods. The accuracies of 
energies and densities when using these methods are also considered. 
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I. INTRODUCTION 



Entanglement, one of the stranger features of quantum me- 
chanics, is now seen as a resource that can be exploited 
for teleportation of quantum states and secure distribution of 
cryptographic keysii In addition, it is thought of as one of 
the main reasons that quantum information devices may be 
able to outperform their classical counterparts. In this respect 
many solid state systems have merit as quantum information 
processors, for example quantum dots.^ However, modelling 
solid state many-body systems exactly is often computation- 
ally intractable, hence approximations are used. One effi- 
cient method to model these systems is to use approximations 
within density-functional theory (DFT).^ DFT maps an inter- 
acting many-body system onto a non-interacting one, from 
which the ground state properties of the former can be cal- 
culated in principle exactly. In practice though, as the key 
quantity in DFT — the exchange-coiTelation potential — is an 
unknown functional of the density, approximations must be 
used. If we consider the entanglement of the ground state, 
then a relationship between DFT and the entanglement must 
exist;^"^ in this paper we propound ways to calculate the 
ground state entanglement of many-body systems using DFT. 
We will consider the entanglement generated by the spatial 
degrees of freedom, and use Hooke's atom as a test system. 
Hooke's atom represents a possible model for describing two 
electrons trapped in a quantum dot, a system similar to the 
ones proposed in various quantum information applications. 
The interacting wave-function of this system can be calculated 
exactly.^ This allows us to quantify its exact entanglement by 
calculating the Von Neumann and linear entropies. We then 
propose two possible ways to calculate the same entangle- 
ment using DFT. This is a delicate issue, as, by construction, 
DFT reproduces the exact density of the many-body system, 
but not its wave-function. To calculate the entanglement, it is 
then necessary to define a suitable 'interacting wave-function' 
within the DFT scheme. As a different route, we also explore 
the possibility of calculating the entanglement of this many- 
body system by means of perturbation theory. We wish to 
discover how well these approximations reproduce the entan- 
glement of the system and therefore ascertain whether they 
are good methods to model a system's suitability as a quan- 
tum information device (e.g. as a component for a quantum 



computer). 

In Section II we introduce the system: Hooke's atom. Sec- 
tion III discusses measures of continuous entanglement and 
uses them to quantify the exact entanglement for our sys- 
tem. We also use the expectation values of the potential and 
Coulomb energies to express the degree of interaction as the 
confining potential varies. In section IV we introduce DFT 
and use the local-density approximation (LDA) to calculate 
the approximate density and examine the link between the ex- 
act exchange-correlation energy and the entanglement. The 
accuracy of the approximate entanglement when using the 
LDA is considered in section V where we define and calcu- 
late approximations to the 'interacting LDA wave-function'. 
In Section VI the second method to link DFT to an interact- 
ing wave-function is proposed. This is based on perturbation 
theory, where the Kohn Sham equations are used as the zeroth- 
order Hamiltonian. We calculate the total energy, the density, 
and the entanglement to first-order in the perturbation using 
both the exact exchange-coiTelation potential and its LDA ap- 
proximation; these results are then compared with first-order 
standard perturbation theory. Finally, the paper concludes in 
Section VII with an overview of the different methods' abil- 
ity to reproduce the exact entanglement and ideas for future 
work. 



II. THE SYSTEM 

The spatial entanglement of interacting electron systems 
has recently been investigated and calculated, by Osenda and 
Serra** using the Von Neumann entropy, for the spherical 
helium model where the Coulomb repulsion is replaced by 
its spherical average. Here we consider a system with full 
electron-electron interaction. To do this we use the model 
system of two interacting electrons in a harmonic potential: 
Hooke's atom. This system may be solved exactly^ and 2D 
and 3D harmonic potentials have been used to model particles 
trapped in quantum dots, for example see Refs.lsil-fTsl 

In atomic units [e = h — m — I/Atteq = 1) the Hamilto- 
nian for two electrons confined by a harmonic potential is 
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where uj is the characteristic frequency of the confining poten- 
tial. 

This system has full electron-electron interaction yet may 
be separated into relative motion and centre of mass co- 
ordinates then solved exactly using the method of Taut7 This 
entails choosing cj so that the three term recurrence relation 
resulting from using a power series expansion for the solution 
of the relative motion radial part terminates. Unfortunately 
this means that exact ground-state solutions exist only for cer- 
tain UJ, the largest of which is 0.5. For other values of uj we 
may construct an approximate solution to the relative motion 
part by using a finite number of terms (400) in the now in- 
finite power series and numerically finding the eigenvalue of 
the relative motion part that causes its wave-function to tend 
to zero at large r. 



III. MEASURES OF ENTANGLEMENT AND EXACT 
RESULTS 

Here we are considering fermionic entanglement. For the 
two electrons the entanglement is spread over the spin and 
spatial parts of the wave-function. We consider the ground 
state of the system so the spin part is that of a singlet and 
thus always maximally entangled, hence we concentrate on 
the spatial degrees of freedom. The entanglement generated 
by the continuous spatial degrees of freedom will depend upon 
the interplay between the mutual electron repulsion and the 
strength of the confining potential. Therefore, in this paper, 
we study how the spatial entanglement changes with the sys- 
tem parameters. We shall quantify the entanglement using 
different measures i.e different entanglement entropies. 



A. Linear entropy 

The linear entropy of the reduced density matrix 

L = Trpi-ed - Trp2 d = 1 - Trp^^^ 

may be used as a measure of entanglement for a pure state (for 
example in Ref. 5) by giving an indication of the number and 
spread of terms in the Schmidt decomposition of the state. 

We quantify the spatial entanglement of the electrons in our 
system by using the analogue of L in the continuous case, 
where 
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FIG. 1: Linear entropy of the reduced density matrix versus u}. 

The numerical calculations for varying lu are displayed in 
Fig- [11 Here we can observe, firstly, that the spatial entangle- 
ment is very low. The linear entropy of the reduced density 
matrix is bounded by d/{d — 1), where d is the dimension 
of the state. As the spatial part is continuous then d is infi- 
nite; hence the bound is unity, which our results are all very 
much less than. The entanglement increases with decreasing 
confining potential as this means that the relative strength of 
the electron-electron interaction increases. According to the 
value of w, we can identify three different regions. For large 
values of the entanglement is essentially zero as we can 
neglect the electron-electron interaction compared to the con- 
fining potential. In essence, in this region, we have two 3D 
isotropic oscillators hence the wave-function is a product in 
electron co-ordinates. Another way to see this is to consider 
that as the confining potential tends to infinity then the elec- 
trons will be forced to overlap. This is allowed for a singlet 
state and the spatial part may be written as |0) |0) which is 
clearly a product state. For intermediate values of lo, we see a 
sharp increase of the entanglement which corresponds to the 
confining potential and electron-electron interaction becom- 
ing of similar magnitude. For lj < 0.001 then the increase of 
entanglement per order of magnitude diminishes significantly. 

We note that a quantum dot of width 2A = lOnm corre- 
sponds to a confining potential of a; = 3.65 Hartrees* for 
GaAs and uj = 0.684 Hartrees* for CdSe. Here Hartrees* 
are 'effective Hartrees' calculated using the related material 

parameters and we have used the formula \ = \ — ^ . The 
entanglement is therefore equivalent to the entanglement of 
Hooke's atom at cj = 3.65 Hartrees for GaAs and lu = 0.684 
Hartrees for CdSe. 



PKd{ri,r2) = y ^'*(ri,r3)*(r2,r3)dr3, 



1. Interaction regions 

To aid understanding of where the different regions of in- 
teraction occur we present a plot of the ratio of the expec- 
tation of the Coulomb interaction to the expectation of the 
external potential (Fig. |2l). This allows us to identify three 
regions that help explain the behaviour of the entanglement in 
Fig. [U for ijj < 0.001 we have that the ratio is greater than 
1.5 and we deem this the 'strongly interacting region' where 
the logarithmic plot of the entanglement appears to plateau; 
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2 > w > 0.001 corresponds to an 'intermediate interaction 
region' where the logarithmic plot of the entanglement in- 
creases sharply; for w > 2 the ratio shows that the external 
potential is the dominant term and we are entering the essen- 
tially non-interacting region where the entanglement is close 
to zero. 
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FIG. 2: Ratio of the expectation of the Coulomb interaction to the 
expectation of the external potential versus u. 



B. Von Neumann entropy 

The Von Neumann entropy of the reduced density matrix is 
given by 

S = -Trpred l0g2 Pred = - ^ A Jo52Aj, 

i 

and quantifies the entanglement for a bipartite pure state. S is 
thought of as the definitive measure of entanglement for a pure 
state (see for example Ref. [3 and the references contained 
therein). Here are the eigenvalues of the reduced density 
matrix. 

To calculate the continuous spatial entanglement, we first 
note that in the discrete case the eigenvalue equation may be 
written as 

^Pred,y0j = A0i. 

j 

Hence in the continuous case we have 

J Pr<:d{iJ)(f>{j)dj = (1) 

This may be solved by turning it into an algebraic problem, 
following the idea of Osenda and Serra^^ 
We use a basis for the eigenfunction 

n' 

and multiply Eq. ([T]i on the left by ?7„ (r) and integrate to give 
the eigenvalue problem 

^(A„„- -\M.nn')Cn' =0. 
n' 



Here 

Ann' = J Vn{r)p'^'^{r,r')r,n'{r')drdr' 

and 

Mnn' = jT]n{r)rin'ir)drdr'. 

To reduce calculation time and increase accuracy an orthonor- 
mal basis is created such that 



Att rii{r)i]j{r)r'^dr ^ 5ij. 

JO 

This is done by using the Gram-Schmidt process; we create 
polynomials Pi{r) so that 

POO 

/ Pi{r)pj{r)w{r)dr = 
Jo 

with the weight function chosen as w{r) = r^e~'^'^^^, 
where ujr — We may then construct rii{r) = 

Pi{r)e~'^''^ 1"^ and then normalise to give the orthonormal ba- 
sis set. The form of r\i (r) is chosen to be similar to the exact 
solution. This process gives the results (Fig. [3]) which are of 
a similar shape as the linear entropy (Fig. [U but of higher nu- 
merical values. This is because for an infinite dimensional 
system, the Von Neumann entropy is unbounded unlike the 
linear entropy. We compare the two measures of entangle- 

0.14^ — ^ — — ^ — ^ — ^ — ^ — ^ — ^ 
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FIG. 3: Von Neumann entropy of the reduced density matrix versus 
u. 

ment (Fig.|4]l by scaling L by 3.75 to achieve equality of the 
measures for the smallest value of w plotted. The measures do 
not scale exactly onto each other, but it is reassuring to see that 
they coincide in behaviour, at least for this system. It appears 
that, for this confining potential range, the easier to calculate 
L, suitably scaled, is a good approximation to S. 

C. Position-space information entropy 

The position-space information entropy Sn has been used 
to study the Moshinsky atomJ^ while the information entropy 
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FIG. 4: Comparison of S and L scaled by 3.75 versus uj. 
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FIG. 5: Position-space information entropy Sn versus oj. 



of the pair correlation function for Hooke's atom has also been 
studied by Atre et al,'* The position-space information en- 
tropy 



Sn = 



n{r) \nn{r)dr 



(2) 



may also be thought of as a zeroth-order approximation to S, 
where the off-diagonal terms of the reduced density matrix are 
set to zero. Let us write 



Pred Pdiag i Poff diag 



where 



n{r) 



n is the electron density and N is the number of electrons. 
Then, when neglecting off-diagonal terms we have 

S ^ - J pdiag log2 Pdmgdr 

= / n{r) mn{r)ar ~\ 

iVln2 / ^ ^ ^ ln2 
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Although much easier to calculate, our results in Fig. |5] show 
that the features of an exact measure of the entanglement are 
lost completely when using the position-space information en- 
tropy as an approximation to S. 



IV. DFT TREATMENT OF HOOKE'S ATOM 

Solving the Schrodinger equation directly for an interact- 
ing many-electron system is computationally onerous. One 
framework that offers a more efficient method is density- 
functional theory, where the density rather than the many- 
body wave-function is employed to calculate properties of the 
system. The theorem of Hohenberg and Kohn'^ states that 
for a non-degenerate ground state of an interacting electron 



system, the density uniquely determines the many-body wave- 
function. Hence in theory all aspects of the ground state of the 
system are derivable from it. The theory predicts a mapping 
between two systems having the same ground state density: 
the interacting system and the 'Kohn Sham' system formed 
by non-interacting particles. The density is found, exactly 
in principle, by solving single-electron equations describing 
the non-interacting system: the Kohn Sham (KS)'" equations. 
They contain the functional derivative of the exchange corre- 
lation energy functional iJ^c which allows us to move between 
the interacting and non-interacting systems. 
The KS equation for our system is 



1 2 2 
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+ VH + v^c{r;[n]))<j){r) ^ e(f>{r). (3) 



Here 



is the exchange-correlation potential, and 

n(r2) 



rdro 



\ri - r2 



(4) 



is the Hartree potential. 



A. Exchange-correlation energy as an entanglement indicator 

Hooke's atom is one of the few systems where the exact E^^ 
may be calculated by inversion of the KS equation.— We use 
this method to calculate the exact Ej^c for Hooke's atom for a 
range of u. The results are depicted in Fig.|6] 

As all the information about interactions is contained in 
i?xc, this would be expected to have a relationship with the en- 
tanglement. It can be seen that i?xc increases towards zero as uj 
decreases but appears to be unbounded as uj increases. Hence, 
it bears little similarity to the entanglement behaviour. This is 
because the interaction energy must be placed in context of the 
total system energy; therefore a better indicator would be E^^ 
as a fraction of the energy of the system (Fig.|2)- This follows 
a similar shape to the entanglement plot (Fig. [Tli: it captures 
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FIG. 6: Exact iS^c versus lj. 



the general trend of the entanglement in the different interac- 
tion regions, but the increase in the 'intermediate interaction 
region' is less sharp and the graph plateaus more slowly in 
the 'strongly interacting region'. Hence, we can think of this 
quantity as an indicator of the degree of entanglement for this 
system. 




mating the KS equation using central differences with the as- 
sumption that the ground state density is spherically symmet- 
ric. In Fig. [8] we present a comparison of ni^DA{r), the self- 
consistent density resulting from using the LDA, and n{r), 
the exact density, for three different values of uj, as labelled. 
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FIG. 7: \E^c/E\ versus as a possible entanglement indicator. 




B. LDA density results 

Unfortunately, the general functional form of _Exc is un- 
known hence approximations are used. One of the simplest, 
yet often effective, is that of i?xc for a homogeneous elec- 
tron gas, this is known as the local-density approximation 
(LDA)."* 

To aid understanding of how the LDA approximates entan- 
glement we first investigate how well the LDA reproduces 
the effects of many-body interactions for Hooke's atom. To 
achieve this we solve the KS equation (Eq. (O) with this ap- 
proximation: for Vxc we use u^^^ = v}f^ + v^;^^ where v^^^ 
is the exchange potential and v]:^'^ is the correlation poten- 
tial. We employ the fit of Perdew and Wang^*^ to Monte Carlo 
simulations for the correlation part while the expression for 
the exchange part is taken from Wigner and Seitz.^' The KS 
equation is solved iteratively until self consistency is reached 
for the density n{r) = 2|(/)(r)p. Our solution is achieved by 
diagonalising the tri-diagonal matrix resulting from approxi- 
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FIG. 8: Comparison of the LDA density (dashed lines) and exact 
density (solid lines) plotted against r for three different values of uj. 

We quantify the relative distance between the LDA and ex- 
act density as 

n % Error = ioo2>i^^M^^. (5) 

This is plotted in Fig. |9] 

As ll! decreases then the relative strength of the electron- 
electron interaction increases (see Fig. |2]i. We find that the 
LDA becomes less accurate as the relative strength of the 
many-body interactions in the system increases and that it 
tends to underestimate the density for small to medium r. This 
is corroborated by the results of Taut, Ernst and Eschrig. — 
This is somewhat counter-intuitive as with smaller lu the den- 
sity is more slowly varying and hence more similar to that 
of a uniform electron gas. It should be noted, however, that 
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FIG. 9: Relative error in the LDA density compared to the exact 
density versus ui 

the self interaction, whereby the mean field approach of the 
Hartree potential Eq. (HI means that an electron is subject to 
an interaction with its own density, is not cancelled in the LDA 
approach. This anomalous self interaction will be pronounced 
for this system as there are only two electrons. The effect of 
the self interaction can be observed by the propensity of the 
LDA density to underestimate the exact density at the origin: 
the electrons feel more repulsion than they should so the den- 
sity is spread further out. 



V. ENTANGLEMENT IN DFT: METHOD I 
(THE 'INTERACTING LDA WAVE-FUNCTION') 



reproduces the DFT density. We will do this for Hooke's atom 
within the LDA and name the resultant wave-function as the 
'interacting LDA wave-function' V'LDA.int- We may then cal- 
culate its entanglement. 

Our assumption is that if the LDA density closely approxi- 
mates the exact density then i/'LDA.int will be close to the exact 
wave-function and consequently the same will be true of the 
entanglement. The way we propose to achieve this mapping is 
to consider the LDA KS equations as exact KS equations for a 
fictitious interacting system, characterised by an external po- 
tential v^^^^. To circumvent the problem^^ of the unknown 
exact form of v^^/^, we choose then to design a method which 
reproduces the LDA density and minimises the functional 

Q= (5-1 f + feel*) (6) 

which appears in the KS theorem. This ensures that the wave- 
function is indeed the ground-state as the method is equivalent 
to Levy's constrained searchi^l Here T is the kinetic energy 
operator and V^e is the Coulomb energy operator. Therefore, 
of the wave-functions which reproduce the density we choose 
the one that minimises Eq. (|6]): this is the wave-function that 
minimises the total energy of the system as the energy arising 
from the external potential, / v^^^'^{r)n]jot^{r)dr, is fixed by 
the density. 

Motivated by the form of the exact solution for Hooke's 
atom we choose the form of the wave-function we use in the 
search to be separable in centre of mass and relative motion 
co-ordinates 



DFT shows that all ground-state properties, therefore the 
entanglement, may be written as a functional of the ground- 
state single-particle density n{r). Unfortunately the explicit 
form of this functional for many properties, including the 
entanglement, is unknown. However we know how to ex- 
press the entanglement as a function of the interacting wave- 
function. In this section we will propose a way of defin- 
ing such a wave-function within the framework of DFT. To 
test it we will calculate the coiTesponding entanglement of 
Hooke's atom and use the LDA to approximate the exchange- 
coiTelation potential. 

We exploit the fact that the density uniquely determines the 
interacting wave-function, and search for the ground state in- 
teracting wave-function that reproduces the LDA density. 

A. Defining and finding an 'interacting LDA wave-function' 

There are properties of the ground state of an interacting 
electron system which may be expressed easily in terms of 
the ground state wave-function but whose expression in terms 
of the density is unknown. Two example of such properties 
are the kinetic energy and the entanglement. The Hohenberg- 
Kohn theorem.'^ showed that the ground state density for an 
interacting electron system uniquely determines its interacting 
wave-function. Hence, knowing the ground-state density from 
DFT calculations, we can invert the problem and use the DFT 
density to find the ground-state interacting wave-function that 



V'LDA,EAi(r,i?) =p(r)e-5"'-'^'e-2-««', (7) 

■i 

Here w^, w^, and the are the parameters that we vary in the 
search and the subscript EA\ signifies evolutionary algorithm 
(see below). The ^lda,eai(?', R) trial form for V-'LDA.int allows 
us to calculate the expectation Eq. ^ quickly. Moreover, this 
method avoids solving the Schrodinger equation each time or 
being overly restrictive in the form of w^^?^. To optimise the 
parameters we combine an evolutionary algorithm with a gra- 
dient descent algorithm to minimise the difference between 
the LDA density and the density nLOA.int arising from the trial 
interacting wave-function 

f Ei \nU)A{ri) - nLDA,mt(ri)| 

/ — 7 — X ■ \°) 

The outline of the first method to approximate the interacting 
wave-function is described in appendix A. We also consider 
the most general form of a two electron ground state wave- 
function arising from a spherically symmetric potential: 

■4^U)A,EAi{ri,r2,0) = ^ aijkrii{ri)rij{r2)Pk{cos{9)), (9) 

ijk 

where 

V^{r) ^ P^{r)e'^'^''" (10) 
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are the functions orthonormal with respect to spherical polar our approximate V'LDA.int to be sufficiently accurate only when 
integration detailed in Section ITlI Bl and we use 64 terms. In the relative error Eq. ^ is less than a few percent, 
this method we consider the sum of / and Q as the function 
to minimise. 

1. Entanglement results 



B. Approximation to the interacting LDA wave-function 
results 

Fig.fTOlcompares nu:,\{r) with nLDA,int(''), the density aris- 
ing from V-'LDA,EAh for the best density match and lowest ex- 
pectation Eq. (|6]l the algorithm found. The error in match 
Eq. (O is displayed also (inset of Fig.fTOli. The match becomes 
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FIG. 10: Comparison of the LDA density (dashed line) to the density 
arising from i/jlda,eai (solid line) plotted against r for three different 
values of uj. 



The entanglement quantified using L for the exact wave- 
function and the approximations to ?/'LDA,int is depicted in 
Fig. [TTl for values of uj that yielded a sufficiently accurate 
approximation to V'LDA.int; i-e / < 4.1% for i/'lda.eai and 
/ < 1% for V'LDA.EAi ■ The next smallest lo for which an 
analytical exact solution to the Hooke's atom problem may 
be found led to / w 8% hence was not included. It can be 
seen that as w < 2 ('intermediate strength' coupling regime) 
then the entanglement of our approximations to f/'LDA.int is 
greater than the true entanglement. We find a large discrep- 
ancy between the entanglement when using the approxima- 
tions to T/iLDA.int and the exact wave-function when the inter- 
action strength becomes comparable to the confining potential 
see Fig.|2] Our calculations show that the entanglement calcu- 
lated from this method is very sensitive to the trial form used 
for the wave-function and to the details of the fitness function 
(/ or / + Q) used. 

For larger uj the entanglement is reproduced. We note that 
for ijj = 1, / is as low as 0.027% and the entanglement is 
still overestimated using i/'lda,eai- However when we con- 
sider '!/'LDA,EA2 for the intermediate regime the match is gener- 
ally better and the overestimate is reduced. Carrying out these 
calculations we realised that entanglement is very sensitive to 
small changes in the wave-function. We have found approx- 
imations to i/'LDA.int that reproduce the tilda accurately, have 
as low a energy as possible and are ground-states by virtue 
of their lack of nodes; however we can not guarantee that the 
energy is indeed the minimum. The change in entanglement 
upon considering a more general basis set and so, in general, 
increasing the accuracy of the match, confirms for example 
that the energy from i/'lda.eai is indeed not the true ground 
state. Due to this high sensitivity of the mapping between den- 
sity and corresponding wave-function and of the entanglement 
on the wave-function details we can not be confident that we 
have found a close enough approximation to the exact V'LDA.im 
to represent the entanglement trend correctly. We speculate 
that a much larger number of iterations of the algorithm could 
improve accuracy, but this would become too computationally 
expensive to make this method a practical way to approximate 
the 'interacting LDA wave-function' . 



C. Position-space information entropy 



less accurate as the importance of many-body interactions in- 
creases. As we based this search around a wave-function of 
similar form to the exact solution then it is not surprising that 
the match is less accurate when the interactions increase as 
this is when the LDA density has a pronounced difference to 
the exact density. It still should be noted that all these matches 
are much better than the match between the exact density and 
LDA density, compare Fig. [8] and Fig.[TO] We shall consider 



We return to the position-space information entropy Sn 
Eq. (|2]l and calculate it using wlda- In Fig. [12] it can be ob- 
served that in this case the LDA gives almost the same infor- 
mation entropy as the exact solution. It is interesting that even 
when n{r) and nLDA(?') are significantly different, i.e when 
ljj is very small, the corresponding LDA and exact position- 
space information entropies are still very similar. Clearly 
this measure is not sensitive to idiosyncrasies in the density 
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FIG. 11: Comparison of L for the exact wave-function and the 
approximate interacting wave-functions (V'lda.eai and '(/'lda,ea2) 
against u. Inset: relative error in the density match / versus oj 
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FIG. 12: Position-space Information Entropy Sn for the LDA density 
(crosses) and the exact density (solid line) plotted against co. 



shape. Sn is a relatively simple integral of a normalised quan- 
tity, the density. As such, if two normalised functions are 
non-negligible in the same region and have roughly the same 
shape — as n{r) and tilda do (see Fig.[8j — they will produce 
very similar values for Sn- This insensitivity is more evidence 
that Sn is a poor approximation to an entanglement measure. 



\ri - r2\ 



(13) 

is treated as the perturbative term. Then the solution to 
Eq. (ITTl l may be treated as a perturbation of the KS equa- 
tions. When Uxc has to be approximated, n{r) is found self- 
consistently for the ground-state and then used in the expres- 
sion for Vxc for the calculation of the excited states. We 
expect this method to reduce the magnitude of the pertur- 
bation as many-body interactions are already partially ac- 
counted for in the zeroth-order term. We therefore expect this 
to produce more accurate results in respect to a more stan- 
dard perturbation expansion. Similar procedures have been 
implemented^^'^*' to give an exact perturbation expansion for 
i?c. In the following we will apply this method up to first or- 
der and examine how well it reproduces the entanglement of 
the system. We will first perform our calculations by using the 
exact Wxc- Hooke's atom is one of the few systems for which 
Vxc may be expressed exactly by inverting the KS equationsJ^ 
This corresponds to a test of the perturbation scheme only, as 
we have eliminated possible additional approximations in Wxc- 
The details of the calculation of the first-order wave- 
function are given in Appendix B. 



A. Comparison with standard perturbation theory 

We may also carry out a standard perturbation on the 3D 
harmonic oscillator equations. Here we define 



where 



and 



H = Ho + H' 



i=l,2 



VI. ENTANGLEMENT IN DFT: METHOD II 
(PERTURBATION USING THE KS EQUATIONS) 

The second method we propose to calculate entanglement 
within DFT, uses the KS equations as a zeroth-order Hamil- 
tonian in a perturbation scheme chosen to give the exact 
Schrodinger equation. This scheme gives another method of 
associating an interacting wave-function with DFT. 

We write the full two-electron Hamiltonian for Hooke's 
atom as 



\ri - r2 



We use a similar method to that described in Appendix B 
to calculate the wave-function to first order. We analytically 
calculate the total energy to first order as 



H = Ha + H' 



(11) 



B. Energy results 



where 

Ho= J2 -l^^ + l^^''^ + N) + in; [n]) 



(12) 



We first compare the results for the total energy of the sys- 
tem. In Fig.[T3]we plot the relative error 



Energy % Error = 100 ( ^^^p" ^ 

E 



(14) 
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where E is the exact system energy and -Eappx is the system en- 
ergy calculated using an approximation: the LDA; first-order 
perturbation of the exact KS equations; first-order perturba- 
tion of the LDA KS equations; or first-order standard pertur- 
bation. It can be seen that standard first-order perturbation 
fares poorly in calculating the energy, while this is greatly im- 
proved by the perturbation of the LDA KS equations to first 
order but only marginally improved again by the perturbation 
of the exact KS equations. The LDA energy though is the 
most accurate and performs very well with an error always 
less than 6% for the plotted confining potentials. It is interest- 
ing that the LDA reproduces the energy much better for weak 
confining potentials than it reproduces the density or entan- 
glement. As an additional check we may calculate the energy 
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FIG. 13: Relative error of the approximate energy from LDA, first- 
order perturbation of the LDA KS equations, first-order perturbation 
of the exact KS equations and standard first-order perturbation, com- 
pared to the exact energy plotted against cj. 

to second order for the perturbation of the LDA KS equations 
(Fig.O. This can be seen to reduce the error, especially for 
the low confining potentials, although it still does not achieve 
as high an accuracy as the LDA energy. In the inset of Fig. [14] 
we show that when implementing a perturbation of the LDA 
KS equations the first-order energy is generally more accu- 
rate than the zeroth-order energy. However, the non-trivial 
behaviour of E''"^ (twice the KS eigenvalue) results in propi- 
tious closeness of E''"' to the exact energy for oj ~ 0.00584. 



C. Entanglement results 

Fig. [15] shows that standard perturbation approximates the 
entanglement well in the essentially non-interacting region but 
poorly for confining potentials smaller than oj w 0.5. The 
perturbation of the KS equations to first order performs much 
better both for the exact v^c and its approximation v]^^^. It 
reproduces the entanglement well for uj > 0.5 after which it 
begins to overestimate the entanglement. This overestimate 
becomes worse as u decreases, but the accuracy is still supe- 
rior to that of standard perturbation. An explanation for the 
behaviour of the approximations to the entanglement shown 
in Fig. [15] is that at strong confining potentials the perturba- 
tion is very small and hence first-order perturbation will give 
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FIG. 14: Comparison of the relative errors of first and second-order 
perturbations of the LDA KS equations together with the relative er- 
ror in the LDA energy plotted against co. Inset: Comparison of the 
relative errors of zeroth (dashed line) and first-order (solid line) per- 
turbations of the LDA KS equations compared to the exact energy 
plotted against oj. 



an accurate wave-function with the correct entanglement both 
for the standard perturbation and perturbation of the KS equa- 
tions. When the interaction increases perturbation theory is 
bound to become less accurate. In particular the resulting first- 
order wave-function will be mainly composed of a mixture of 
excited state wave-functions, therefore the entanglement will 
increase. We find that this mixture consistently overestimates 
the entanglement. 

Our hypothesis was that the perturbation term H' 
(Eq. fT3[). when using the KS equations as the zeroth-order 
Hamiltonian, was smaller than the standard perturbation term, 
due to some of the interaction being taken into account in the 
KS equations, which represent, in this case, our zeroth-order 
Hamiltonian. This implies that the related entanglement will 
be more accurate, which is supported by our results. To in- 
vestigate the accuracy of the first-order perturbation of the KS 
equations when v^c is approximated we also consider this per- 
turbation using v}^^^. We see that the entanglement calculated 
using perturbation of the LDA KS equations is very similar 
to that calculated with perturbation of the exact KS equations, 
although there is a slight deterioration in the accuracy. There- 
fore the LDA may be thought of as sufficiently accurate for 
this method and little is gained by using the exact exchange- 
correlation potential. Hence for smaller confining potentials 
(cj < 0.5) it seems that higher-order perturbation is required 
to improve the accuracy of the calculated entanglement, not a 
more accurate exchange-correlation potential, at least for this 
system. Around ut « 0.01 the perturbed LDA KS entangle- 
ment starts to decrease slightly; after this it appears that the 
numerical accuracy becomes insufficient. The range we can 
accurately perform the calculation on is similar to the range 
we were successful in approximating V'LDA.int, hence it allows 
us to compare the two methods. 

We may also quantify the accuracy of the method by com- 
paring the match between the density calculated using the var- 
ious approximations and the exact density, using an expres- 
sion similar to Eq. ([5]). As can be seen in Fig.[T6]the perturba- 
tion using the exact KS performs slightly better than uloa so 
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FIG. 15: Entanglement quantified with L of the exact solution, first- 
order perturbation of the KS equations, first-order perturbation of 
the LDA KS equations and standard first-order perturbation, plotted 
against uj. 



reversing the result of the energies. The standard perturbation 
is poor and becomes rapidly worse as lu decreases. We also 
compare the perturbation of the exact KS equations to the per- 
turbation of the LDA KS regarding the accuracy of the density 
in the inset of Fig. [16] As with the entanglement there is an 
improvement when using the exact Wxc, especially for small lu, 
but the increase in accuracy is much larger for the density. 
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FIG. 16: Comparison of the relative error in the reproduction of the 
density by approximate methods versus u. Inset: comparison of the 
relative error in the reproduction of the density between perturbation 
of the exact KS equations and perturbation of the LDA KS equations 
versus u). 



VII. OVERVIEW AND CONCLUSIONS 

We have calculated the entanglement for the spatial degrees 
of freedom of two electrons confined by a harmonic potential 
(Hooke's atom) and compared different measures of the en- 
tanglement for this system. We found that the linear and Von 
Neumann entropy of the reduced density matrix give consis- 
tently similar results: for Hooke's atom in its ground state 
the entanglement increases with increasing electron-electron 
interaction but is always quite low. The behaviour of the en- 
tanglement, in this case, is therefore independent of which of 
these exact measures is employed. We have also studied the 



position-space information entropy, which does not produce 
the coiTect behaviour of entanglement when the strength of the 
external potential is varied. As the position-space information 
entropy can be viewed as a simple and appealing approxima- 
tion to the Von Neumann entropy we conclude that caution 
must be used when dealing with it; more analysis is neces- 
sary to understand for which systems such an approximation 
is valid. We analysed a possible entanglement indicator, the 
exchange-correlation energy E^^. We found that despite all 
the many-body interactions being contained within this quan- 
tity it does not reproduce the features of the entanglement. We 
calculated \E^c/ E\ and found it a better indicator of entangle- 
ment: it gave a similar trend to that of S or L with respect to 
the confining potential. 

In the second part of the paper we introduced two possible 
ways of calculating the entanglement of a many-body system 
within density-functional theory. This is a powerful method to 
treat many-body systems, but there is no obvious way within 
it to calculate the system wave-function, which we needed 
to deduce the entanglement. To do so we have proposed an 
'interacting LDA wave-function' ^/'LDA.int and a perturbation 
scheme based on the KS equations. We have applied these 
methods to the Hooke's atom problem, the first using the LDA 
to approximate the exchange-correlation potential v^c and the 
second using the exact Uxc and its LDA approximation. 

As a precursor to understanding the entanglement we 
looked at the accuracy of the density reproduction when us- 
ing the LDA. We found that when using the LDA to model 
this system the resultant density is very close to the exact den- 
sity when the confining potential is large, however the match 
becomes progressively worse as the relative strength of the 
electron-electron interaction increases. 

To reproduce the entanglement, we defined an 'interact- 
ing LDA wave-function' V'LDA.int, which reproduces the LDA 
density and minimises the expectation of the kinetic and 
Coulomb energies. As an approximation to 'i/'LDA.int, V'lda.eai 
and i/'LDA,EA2 were found using a combined evolutionary and 
steepest descent algorithm. These approximations were used 
to calculate the entanglement and revealed that they overesti- 
mate the entanglement when there is non-negligible electron- 
electron interaction {uj < 2) (Fig. [TTjand Fig. [TSTl. However 
we found a discrepancy between the entanglement calculated 
using different trial functions even for very accurate density 
matching. This suggests that the mapping between the den- 
sity and the interacting wave-function in a DFT scheme is 
highly non-linear and confirms the sensitivity of the entan- 
glement to the details of the wave-function. As a result we 
cannot know whether the ground-state we have found is close 
enough to that of the interacting LDA system, so this method 
to find the 'interacting LDA wave-function' does not work 
well. We speculate that the search for the interacting LDA 
wave-function would be more successful if it were possible to 
search generally and efficiently for the external potential of an 
interacting electron system that gives the LDA density. 

We then considered perturbative methods. We first analysed 
the behaviour of the system total energy. The results showed 
that standard first-order perturbation is inaccurate but is much 
improved by using the exact KS equations as a zeroth-order 
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Hamiltonian, with only a small decrease in accuracy when the 
LDA KS equations are used instead. The LDA energies are 
the most accurate. This shows how errors are cancelled in the 
LDA to give accurate energies despite the density and entan- 
glement being inaccurate for strong interactions. For the ac- 
curacy of the density reproduction we observed that standard 
first-order perturbation performed worst and a large improve- 
ment was given by perturbing the LDA KS equations. This 
gave an even higher accuracy than the LDA density. Density 
reproduction was further improved (markedly so for small cj) 
by using the KS equations with the exact Vxc. 

The results obtained for the entanglement are summarised 
in Fig. [17] and Fig. [18] The entanglement calculated using 
first-order perturbation of the LDA KS equations was larger 
than the exact entanglement but not until after w « 0.5 so this 
is already an improvement on the results from the approxi- 
mations to V'LDA.int- In addition, when it failed to accurately 
approximate the entanglement, the overestimate was less se- 
vere than when using V'lda.eai or V-'lda,ea2- The accuracy of 
the entanglement when using perturbation of the KS equations 
was only slightly improved by using the exact v^c- This sug- 
gested that higher-order terms are required, not more accu- 
rate exchange-correlation energy functionals, when employ- 
ing this method. Standard perturbation was again the least 
accurate amongst the perturbative methods, severely overes- 
timating the entanglement when the confining potential was 
small. This supported the hypothesis that a perturbation of the 
KS equations would approximate the exact results more accu- 
rately than standard perturbation as some of the interaction is 
already included in the KS equations. Even though far from 
optimal, standard perturbation is stiU an improvement over the 
results for the approximations to V'LDA.int for to > 0.03. We 
also note that these trends are observed whether we use the 
linear entropy of the reduced density matrix L or the Von Neu- 
mann entropy of the reduced density matrix 5*, as can be seen 
by comparing Fig. [17] and Fig. [18] 
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FIG. 18: Overview of S for the approximate methods and exact so- 
lution plotted against iu. 



is less clear, though, why the entanglement is overestimated, 
rather than underestimated, when using perturbation of the ex- 
act KS equations and standard perturbation. 

The results have indicated that neither our approximation 
to ^LDA.int nor first-ordcr perturbation approximate the spatial 
entanglement accurately when considering strongly interact- 
ing systems, yet first-order perturbation give the energy to a 
relatively high accuracy. Hence it seems that entanglement is 
much more susceptible to small changes in the density than 
the energy: the entanglement is an integrated quantity of a 
wave-function, which is very sensitive to changes in the den- 
sity. 

The issue of calculating an entanglement within the LDA 
has also been explored by V. V. Franfa and K. Capelle,*" but in 
the context of the Hubbard model and spin entanglement. We 
concur with their conclusion that spatial inhomogeneity (in 
our case determined by the strength of the confining potential) 
has a major effect upon entanglement. 

Further work is planned to investigate whether we can ef- 
ficiently find the external potential of the interacting system 
that reproduces the LDA density: the ground-state of such 
a system would be the true 'LDA interacting-wave-function' 
^LDA.int- We then intend to investigate whether applying a 
self-interaction correction (SIC) to the LDA will improve the 
approximation of the density: the self interaction becomes a 
predominant error as the number of electrons decrease, hence 
by quelling this anomalous effect to some degree in Hooke's 
atom we would hope to make improvements in density accu- 
racy and therefore the approximation of the entanglement. 



FIG. 17: Overview of L for the approximate methods and exact so- 
lution plotted against u. 

We find that all our methods to approximate L and S over- 
estimate the entanglement. For the approximations to the en- 
tanglement involving the LDA, we suggest that, in addition to 
the discussed limitations of the methods, it is also an effect 
of the self-interaction (see discussion of Fig.|8]l, which effec- 
tively enhances the interaction and hence the entanglement. It 
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Appendix A: Combined evolutionary and gradient descent 
algorithm 

The following is the layout of the computational method we 
used to calculate i/iLDA.int- 

• Create a collection of wave-functions with random pa- 
rameters (polynomial coefficients and ujr, i^r)- 

• Minimise the fitness function (Eq. (O) by using a gradi- 
ent descent algorithm on the parameters of each wave- 
function 

• Keep the closest matches and create other wave- 
functions by randomly mutating parameters in the clos- 
est matches and creating new wave-functions with ran- 
dom parameters. 

• The previous two steps are then repeated. 

• After around thirty generations the program ends. Of 
the overall best ten matches, those with the lowest ex- 
pectations of the kinetic and Coulomb energy (Eq. (|6]l) 
are chosen as the best approximations to the interacting 
ground-state wave-function. 

To ensure that we have found a ground-state wave-function 
the ten best matches are also checked to verify that they are 
node-less and that they are lower in energy than the non- 
interacting wave-function formed by the product of the single 
electron KS orbitals. 



Appendix B: First-order perturbation using the KS equations 
We require the ground state wave-function to first-order 

,(1) 



To find ^'n ' we use the Rayleigh-Schrodinger method, where 
the first-order wave-function is written in terms of the unper- 
turbed basis 



We write a single electron wave-function as 

(pnl7n{r) = Rnl{r)Yim{0,(l)). 

Using the ansatz i?„/ (r) = Uni (r) / r gives the radial equa- 
tion 



/ 1 ^2 1 

H — — ]Unl(r) = enlUnl{r). 



2r2 

For I = this is solved self-consistently when t;xc has to be 
approximated, then the ground-state self-consistent density is 
used in the calculation of v^c for the higher angular momen- 
tum states. From the eigenvalues there appear to be no degen- 
eracies due to the non- trivial nature of v^c and vh- 

The exact wave-function is a singlet, therefore any triplet 
states will have zero contribution or will cancel out; hence in 
our scheme we chose to neglect their contribution in Eq. ( IBll ). 
The spin singlet excited state wave-functions can be divided 
into two forms: 

• Product states 



(t)nihmAri)(t)n2l2m2ir2); 

Non-product states 

H^n\l\m\ Vl )0n2l2^2 v^2 ) ~r (Pn2l2^2 v 2 )4'n\l\rai (n)). 

(B3) 



V2 



For the wave-function to contribute to the first-order pertur- 
bation calculation then 



(0) 



H' 



We shall use this to determine the quantum numbers we need 
to consider and hence simplify the numerical calculation. 



Here 



E 



(0) 



k=0 



(0) \ 



H' 



(0) 



(Bl) 



(B2) 



and Qnn — (normalisation of the first-order wave- 
function). Here we concern ourselves with the ground state 
(n = 0) only. We solve 

by separating Hq for each electron, then separating into ra- 
dial and angular parts and solving the radial part using exact 
diagonalisation techniques. 



1. Simplifying the integrals 

a. Product states 

We first consider the fxc contribution from the product 
states (j>nihmi {ri)(f>n2i2ni2 (^2)- Then w^c enters into the per- 
turbation calculation as 



fooo 



(ri).^ooo(r-2)| V ^;.?^(rO !</> 



1=1.2 



rlrldridr2 / Fgo^imirff^i / yoo^i2m2'^^2- 



13 



From the orthogonality of the spherical harmonics then 
mi = TO2 = h = h = for this to be a non-zero integral 
and from the orthogonality of the ground state radial solutions 
ni = 712 = as well. Hence there are no contributions from 
the term for excited states of this form. This result also 
applies to the term arising from vh- 

Next we consider the Coulomb term, which may be ex- 
panded as a series of spherical harmonics: 



where is the angle between ri and r2, we obtain the contri- 
bution due to product states of energy as 

which is symmetric on exchange of positions as required. 



oo +1' 

E E 



An (r<)' 



(B4) 



where r< = min{ri, and r> — max{ri, r2}. We then 



use 



(B5) 



in Eq. iB4i to employ the orthogonality of the spherical har- 
monics of the excited states to the series. This, together 
with the required symmetry on exchange of particles, forc- 
ing Til = 712 and li = I2, means the wave-function must be of 
the form (ri)0„/_™(r2). 

We have, again courtesy of orthogonality, that 



('/>ooo(ri)(?!'ooo(^'2)| 



1 



\ri - r2 



7il—rn 



47r^ ' {21 + 1) 



Rm{r,)Rl^{r2)-^^ 



1 A-jr 

Rni{ri)Rniir2)rlrldndr2 = ^(-1)" (^7:^^) 

where A is understood to be the previous double integral and 
we have used Eq. ( |B5t plus the explicit form of Foo- 

Now as the energy eigenvalues have no m dependence, we 
are at liberty to sum over m when calculating the contribution 
of an excited state. Hence the contribution to Eq. ( IB4b from 
the product states of energy Eni will be 



'm—-\-l 



in 



^(0) _ ^(0) (2Z + 1) An 



{-l)''^Rni{ri)Rni{r2)YU0u<l>i)Yi_,„ie2,cj)2) 



An A 



m—-\-l 



-^(0)_^(0)(2/+T)4^E(-ir^-(^0i?.M 

^0 ^nl ^ ' m=-l 

y,,„(0i,0i)i^:„(^2,02). 

Here we have used Eq. ( IB5l l. Then using that 



ra—-\-l 
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b. Non-product states 

Let us consider for the non-product state (Eq lB3l l the con- 
tribution from the v^^ term, when one electron is in the ground 
state. 

Let us define V'G('^i,'f2) = '/'ooo(»^i)0ooo(^2) and 
ilJE{ri,r2) = {(f>oQo{ri)(t)nini{r2) + (l>ooo{r2)(f>nini{ri)) 
then 



(V'G(n,T2)| ^ Wxc(r'i) |V'is(n,T2) 

1=1,2 

= 2 / / tpGVxciri)'>pEdridr2, 



where we have used the symmetry of the potential on ex- 
changing ri and r2- We note that one of the single electron 
wave-functions must be in the ground state for these integrals 
to be non-zero. By orthogonality the previous equation be- 
comes 

("^gI E "xclj^i) li's) ^ J J '/'OOo(T'l)0OOo(r2) 

nlm {ri)dridr2 

2 



V2 



Rm{ri)vxc{ri)R7ii{ri)r^dri / YQ^YimdU.i, 



m— — / 



where we have used the fact that the wave-function is nor- 
malised. Hence I — m — for this contribution to be non- 
zero, with a similar result for vu- 

For the Coulomb term of this non-product state, the use of 
the expansion in spherical harmonics Eq. (IB4I) forces I = m = 
to ensure that the integral is non zero, with result 

(VgI I ^ 1 li'E) - 4^47r / / i?oo(ri)i?oo(r-2) — 

\ri - r2| \/2 J J 7-> 

{Rm{ri)Rno{r2) + Rm{r2)Rno{ri)) r\rldridr2 — . 

47r 

Next we consider the non-product state where neither electron 
is in the ground state. As neither single particle wave-function 
now has n = and / = 0, then this means that there is no 
contribution from terms with v-^c or vn for this form. For there 
to be a Coulomb contribution, we have that from before li = 
I2 and mi ~ — m2, so we consider only states of the following 
form 

V'S = (0niim(ri)0n2i-m(^2) + f/)™! im (^2 j^nji-m (''l ) ) ■ 
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A similar analysis to that employed in the previous sub- Using the orthogonality of the Legendre polynomials 
section allows the perturbation contribution from wave- 
functions with energy + En^i to be written in terms of 
the radial part and Legendre polynomials: 



B 



^0 ^rnl ^1121 



where 



B 



R*oo{ri)Roo{r2)^^^ 



{Rmi{ri)Rn2i{r2) + Rmi{r2)Rn2i{ri)) r^r^dridr 



Pi{cos{9))Pk{cos{0))dn = 



in 
21 + 1 



Hk 



gives the easily computable expression for the density 



2. Calculating the density and reduced density matrix 

By grouping the wave-functions by we can write to first- 
order 



71 = 2 / /2dri47r + 2 / g^^ny + 2 / e^driy 



*o = /(ri,r2)Po(cos((?)) +5(ri,r2)Pi(cos(0)) + 
h{ri,r2)P2{cos{e)) + ■ ■ ■ 



A similar treatment gives the reduced density matrix. 
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